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ABSTRACT:  Results  are  presented  of  a  computer  simulation  of  the  flow 
induced  in  a  Plexiglass  (PMMA)  rod  by  a  charge  of  tetryl.  The  charge 
and  the  rod  are  the  standard  donor  and  attenuator  of  the  Naval 
Ordnance  Laboratory  Large  Scale  Gap  Test  ( LSGT ) .  The  computations 
provide  a  better  understanding  of  the  flow  and  hence  of  phenomena 
which  have  been  observed  in  PMMA  rods  during  calibration  studies  of 
the  LSGT.  They  also  give  data  on  the  shape  and  duration  of  the 
pressure  pulse  (in  the  PMMA)  which  is  transmitted  into  the  acceptor 
explosive  in  the  LSGT.  Such  data  are  needed  in  the  study  of  shock 
sensitivity  of  explosives. 
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SECTION  1 
INTRODUCTION 


The  NOL  gap  tests  have  long  been  employed  for  characterizing  "high"  or  detonating 
explosives.  In  such  tests  one  usitally  observes  the  limiting  gap  for  detonation  in 
a  rece-tor  by  th>.  shock  wave  transmitted  from  a  donor  charge  through  a  condensed 
inert  medium.  Due  to  the  complex  shock  wave  motion  and  the  short  duration  of 
the  gap  test,  a  complete  understanding  of  the  test  is  difficult  to  obtain  experi¬ 
mentally.  For  a  more  detailed  understanding  of  the  phenomena  in  this  type  of 
test  a  numerical  simulation  of  the  test  should  be  of  value.  report,  prepared 

under  Contract  N60921-69-M-2003  for  the  Naval  Ordnance  Laboratory  and  partially 
funded  under  Honeywell  corporate  funds,  summ?T*izps  the  results  of  the  numerically 
calculated  flow  in  the  simulation  of  the  NOL  regular  card  gap  test. 

The  regular  card  gap  test  consists  of  a  2.  O-inch-diameter  cylinder  of  Plexiglas 

(PMMA)  which  is  shocked  by  a  charge  of  tet^yl  2.  0  inches  in  diameter  and  2.  0 

inches  long.  It  is  described  in  References  1  and  2.  In  the  course  oi  calibrating 

the  test,  several  interesting  phenomena  ar  observed.  One  of  tnese  is  the  "bump'1 

in  the  plot  of  shock  velocity,  U  ,  versus  the  distance  from  the  interface.  This 

s 

phenomenon  is  observed  between  20  and  40  mm  from  the  interface  between  the 
explosive  and  the  PMMA.  Another  interesting  feature  of  the  flow  in  the  PMMA 
cylinder  is  a  discontinuity  in  the  curvature  of  the  snock  front  which  is  observed 
at  about  25  mm  from  the  interface.  Still  another  interesting  feature  of  the  flow 
is  a  sharp  decrease  in  magnitude  of  the  slope  of  the  pressure  at  the  shock  front 
versus  distance  from  the  interface. 

The  purpose  of  the  project  reported  here  was  to  dynamically  simulate  the  explosive 
loading  of  the  PMMA  and  give,  at  least  qualitative,  theoretical  understanding  of 
the  various  phcnor  .ona  observed  experimentally.  This  report  summarizes  this 
simulation  conducted  with  a  :  qrangian  strength  of  materials  code  called  HEMP 
(Reference  3). 
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The  results  presented  here  do  indicate  u  qualitative  understanding  of  the 
phenomena  and  are  thus  encouraging  considering  the  size  of  the  effort  involved. 
However,  it  is  felt  that  they  tall  short  in  several  respects.  First,  the  study 
shows  that  if  it  is  desired  to  obtain  detailed  quantitative  information  of  ffects 
such  as  abrupt  changes  in  radius  of  curvature  of  the  shock  front,  it  will  be 
necessary  to  zone  the  calculations  much  finer. 

Secondly,  if  the  simulation  is  to  yield  detailed  aspects  of  inflections  in  flow 
quantity  graphs  and  the  like,  more  comprehensive  equation-of-state  data  should 
be  obtained.  In  this  regard  it  is  believed  that  more  detailed  information  on 
failure  must  be  obtained  before  that  complicated  phenomenon  can  be  completely 
under'^od.  In  the  last  analysis,  it  must  be  said  that  there  remains  the  possi¬ 
bility  that  a  more  time-consuming  search  of  the  voluminous  data  produced  by 
this  simulation  could  possibly  yield  new  discoveries  unnoticed  by  these  authors. 
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SECTION  2 
THE  HEMP  CODE 


The  HEMP  code  is  a  Lagrangian  hydrodynamic-elasto-plastic  code  developed 
at  the  Lawrence  Radiation  Laboratory  by  Dr.  Mark  Wilkins  and  discussed  in 
Reference  3.  The  basic  equations  of  motion  are  presented  in  what  follows: 

Equations  of  motion  are  in  x-y  coordinates  with  cylindrical  symmetry  about  trie 
x  axis.  (It  is  desirable  to  have  the  problem  formulated  also  in  plane  x-y  coor¬ 
dinates;  for  this  case  the  terms  marked  j  are  set  =  0.  ) 


d  2  d'F  T 

-JS*  +  _JEZ  + 

Sx  hy  y 


3T 


=  Px, 


"Sx 


15^  +  hcL 


*ee  , 


ay 


i 


py , 


(i) 


=  sxx  ■  <p  ^  ■  syy  -  (p  Lee  =  see  '  (p 


(2) 


Equations  of  continuity: 


V.  =  ax  +  3y 
V  *5y 


Energy  equation: 


E  =  -  (p  -  q  >  V  +  V  ( s  e  1  s  e 

1  xx  xx  w  w 


Artificial  viscosity  (quadratic  "a"): 


•>  n  •  > 

q  -  C  “  ;  (V/VP  A/V 


0  0 
where  C  =  con.-  ant.  A  *  zone  area,  and  z 


T  c  ) 

XV  XV 


reference  densitv 


(31 


(4) 


(51 
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Equation  of  state: 


stress  components 


rxy  "  ^£xy*  +  ^xy* 

where  p  =  shear  modulus  and  6  =  correction  for  rotation. 


*xx  ■ 

•„  ’  2H(V-‘/3V/v)+  6yy. 

*ee  -  2^-1/iv/v). 


velocity  strains 


•  _  ax  ;  _  v  ,* 

exx  "  1ST*  e00~T  1 


i  =  e  =  ii+  42 

eyy  ay  xy  ax  ay 


Hydrostatic  pressure: 


p  =  a(p-l)  +  b(p-l)  +  c(t|  D^  +  dpE, 


=  1/V  p/pU 


Von  Mises  yield  condition: 


,  2  ^  2  2.  0/„  ,vo/  „ 

(s 1  +  s2  +  s3  ;  -  2/3  (Y  )  s 


where  Y°  =  material  strength  and  (sr  s2>  Sg)  are  the  principal  stress  deviators. 
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Notation: 


x,  y  soace  coordinates  e  ,  e  .  e  strains 

J  ■  xx  yy  96  xy 

x  velocity  in  x  direction  p  hydrostatic  pressure 

y  velocity  in  y  direction  V  relative  volume 


^  Zyy-  2^0  total  stresses  E 

T  shear  stress  p 

xy 

s  ,  s  ,  s„„  stress  deviators 
xx  yy'  69 


internal  energy  per 
original  volume 

density 


The  dot  over  a  parameter  signifies  a  time  derivative  along  the  particle  path. 


The  finite  difference  equations  solved  in  the  HEMP  code  are  based  on  the 
integral  equations 


3F 

ay 


JcF(n  ‘  j)  dS 


(10) 


where 

C  =  boundary  of  area  A 

S  =  arc  length 

A 

n  =  normal  vector 

A 

t  =  tangent  .  „ctor 


As  such  they  yield  more  accurate  expressions,  especially  in  cases  of  high  grid 
distortion. 
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SECTION  3 

EQUATIONS  OF  STATE,  ZONING  AND  INPUT  DATA 


CONSTITUTIVE  RELATIONSHIP  FOR  PMMA 

Hugoniot  data  listed  in  the  contract  and  discussed  in  Reference  2  gives  the 
shock  velocity  versus  particle  velocity  as  a  piecewise  linear  function 

Us=So  +  SlUp  (11) 

with 


/  S  \ 

r 

1 0.257' 

° 

=  i 

\  Si 

i 

1.61 

\  ll 

i 

for  U  *0.  05 
P 


/  0.295 
1  0.  85 


forU  <0.05! 
P 


with  velocity  measured  in  centimeters  per  microsecond, 
the  Rankin-Hugoniot  equations. 


This  yields,  via 


(12) 


=  1.  185  gm/cm' 


m  =  p/pQ  -  1 


(13) 


lSo\ 

fo.257^ 

1 C.  295^  1 

(sl) 

1.  61 

for  p  >-0.  02, 

^0.85  J  forP  <  °-02J 

(14) 


wherep^  ( p)  is  the  Hugoniot  pressure  and  pq  is  the  initial  density  of 
PMMA.  Pressure  is  measured  in  megabars. 
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The  piecewise  linearity  of  the  shock  velocity-particle  velocity  plot  yields  a 
discontinuity  in  the  derivative  of  pressure  with  respect  to  density  along  the 
Hugoniot.  This  discontinuity  shows  up  on  isentropee  as  well,  and  thus 
affects  the  sound  speed.  Such  discontinuities  in  sound  speed  cause  separations 
in  rarefaction  wares.  The  pressure  Hugoniot  curve  is  shown  in  Figure  1. 

It  is  well  known  that  the  isentropes  for  plastics  diverge  rapidly  from  the 
Hugoniot.  Thus  some  energy  dependence  is  required  in  the  equation  of  state, 
t  or  our  simulations  the  general  form  of  the  Mie-Gru neisen  equation. 


p  -  pH(p)  =  ^  [e  -  eh  (p)] 


(15) 


was  chosen.  Here  y  (p)  is  the  Gruneisen  ratio,  E  is  the  specific  internal 

energy  in  megabars  -  cm  /gm,  t  is  the  specific  volume,  p  is  the  pressure 

in  megabars  for  states  not  necessarily  on  the  Hugoniot,  and  E„(p)  is  the 

n 

Hugoniot  specific  energy 


Er(p) 


!  Ph'®1 

2  P 

o 


m 

1  +m 


(16) 


From  the  viewpoint  of  material  physics  it  is  most  appropriate  to  view 
Equation  (15)  as  a  first-order  expansion  of  the  pressure  about  the  Hugoniot. 
Thus 


3p 


PH 


(P) 


r  (p) 

T 


(1") 


With  this  interpretation,  y(p)  does  not  necessarily  have  the  statistical  mech¬ 
anical  identity  of  the  analogous  quantity  in  Gruneisen's  original  work  on  metals. 
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Since  the  appropriate  value  of  y  for  PMMA  was  not  available  from  NOL 
sources,  data  for  Lucite  from  Wilkins  was  used.  This  value  is 

y  =  0.6  (18) 

With  these  values  inserted  in  Equation  (15)  the  equation  of  state  becomes 

p  =  A(p)  +  B(p)e  ,  e  =  pQE  (19) 

with 

B(p)  =  (20) 

po 

A(p)  =  pH(p)  [l  -  1/2  (0.6)  (p/pQ  -  1)]  (21) 

The  material's  strength  in  tension,  Yq  ,was  taken  from  the  data  of  Keller 
and  Trulio  (Reference  5): 

Y  »  1  kb  (22) 

o 

The  material  was  allowed  to  flow  plastically  after  yield.  The  shear 
modules  of  elasticity,  H-,  was  taken  from  the  Handbook  of  Physics  and 
Chemistry, 

H  =  0. 0143  (23) 


EQUATIONS  OF  STATE  FOR  TETRYL 

A  gamma-law  equation  of  state  was  used  for  tetryl  with  the  detonation  velocity, 
D,  and  the  initial  density,  pQ,  given  by 
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While  this  system  gives  the  correct  detonation  velocity,  initial  density,  and 
energy  release,  it  does  not  neccessarily  yield  the  correct  Chapman-Jouguet 
pressure  and  density.  It  is  however  self-consistent  and  satisfies  the  Rankin- 
Hugoniot  equations  and  the  Chapman-Jouguet  hypothesis.  The  self-consistent 
values  of  the  Chapman-Jouguet  density  p^j  ,  and  pressure  p^  are  then  given 

by 


11 
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-?Ei.  *  .  ,  p~T  =  2.0841068 


a,  D 

PCJ  =  “T  =  0.2156327186 
y  +  1 


Finally,  the  p,  t,  E  equation  of  state  is  given  by 


This  equation  does  not  result  from  the  previous  assumptions.  It  requires 
in  addition  the  assumption  that  all  isentropes  from  the  detonation  Hugoniot 
are  polytropic  with  the  same  polytropic  index  y(i.  e. ,  of  the  form  of  Equation 
(29)  with  the  appropriate  scale  constant). 

It  is  the  opinion  of  these  authors  that  better  detonation  data  could  be  obtained 
if  the  measured  value  of  the  Chapman -Jouguet  sound  speed  cculd  be  obtained. 
This  together  with  the  constants  pQ  and  D  would  then  deteri  line  the  Chapman- 
Jouguet  state.  If,  in  addition,  purely  mechanical  data  on  the  expansion  of  the 
detonation  products  were  available  -  for  example  data  like  that  u'  Lee,  et  al. 
(Reference  7)  -  then  the  detonation  and  detonation  products  would  he  determined 
by  purely  mechanical  means.  At  this  time  no  such  data  exists. 


Finally  it  should  be  stated  that  the  detonation  is  simulated  by  an  energy  dump 
method  (see  Reference  3).  This  method  gives  the  correct  detonation  speed 
but  fails  to  give  a  detonation  pressure  as  high  as  that  given  by  Equation  (31). 
This  point  will  be  discussed  later  in  this  report. 


12 


NOLTR  69-219 


PROBLEM  ZONING 

The  configuration  simulated  (as  shown  below)  was  a  4-inch-long,  2 -inch -diameter 
cylinder  of  PMMA  shocked  by  a  2 -inch-long,  2 -inch-diameter  cylinder  of  ttetryl 
which  was  initiated  centrally  at  the  end.  The  computation  was  carried  out  using 
14  zones  per  inch  (both  radially  and  longitudinally).  This  resulted  in  1, 176 
zones.  The  data  storage  and  computer  program  together  filled  the  core  storage 
of  the  Control  Data  Corporation  3600  computer  with  131,  000  words  of  memory. 
This  computation  was  carried  out  to  a  simulation  time  of  32.  5  microseconds  which 
corresponded  to  about  1.  5  hours  of  computer  time. 


Detonation  Center 
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SECTION  4 

COMPARISON  OF  SIMULATION 
WITH  EXPERIMENTAL  RESULTS 


Comparison  of  the  results  presented  here  with  the  experimental  results  or 
results  obtained  elsewhere  requires  an  explanation  of  how  the  graphs  were 
obtained.  For  computational  purposes,  the  zero  of  axial  distance  was  taken 
to  be  the  beginning  of  the  tetryl.  Thus,  the  interface  between  the  tetryl 
and  the  Plexiglas  was  at  5.08  cm.  Results,  which  are  presented  in  such 
reports  as  NOLTR-65-43,  place  the  zero  of  axial  distance  at  the  interface 
of  tetryl  and  Plexiglas.  Therefore,  a  distance  of  5.  08  cm  had  to  be  sub¬ 
tracted  from  axial  distances  presented  here  in  order  to  compare  them 
with  axial  distances  presented  in  NOLTR-65-43.  Also,  zero  time  in  reports 
such  as  NOLTR-65-43  is  taken  as  the  time  when  the  shock  enters  the 
Plexiglas.  For  the  results  presented  here  zero  time  was  taken  as  the  time 
when  the  detonation  was  initiated  in  the  tetryl.  Determination  of  the  time 
when  the  slock  enters  the  Plexiglas  in  the  numerical  computation  is  resolved 
only  to  within  the  time  step  of  the  computation.  From  the  numerical  results, 
the  shock  entered  the  Plexiglas  in  the  time  interval  of  7.  0  psec  to  8.  0  psec. 
Rather  arbitrarily,  a  zero  time  for  the  shock  to  enter  the  Plexiglas  was 
taken  as  7.25  psec.  This  was  done  so  the  numerically  obtained  shock  loca¬ 
tion  versus  time  curve  initially  coincides  with  that  given  in  NOLTR-65-43, 
Table  A2.  An  explanation  of  now  each  graph  was  obtained  follows. 
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SHOCK  LOCATION  VERSUS  TIME 

The  numerical  computation  was  carried  out  using  the  Q-method  (Reference  3). 
For  this  graph  (Figure  2),  the  shock  location  was  defined  as  the  position  of 
maximum  The  zero  of  the  numerical  results  has  been  adjusted  so  a 
comparison  with  results  from  NOLTR-65-43  can  be  made.  The  zero  adjust¬ 
ment  of  the  numerical  results  has  been  made  by  subtracting  5.  08  cm  from 
each  axial  distance  and  7.25  nsec  from  each  time. 


TIUEpMC 

Figure  2.  Shock  Location  versus  Time 


This  method,  which  seems  to  be  part  of  the  folklore  of  numerical  hydro¬ 
dynamics,  gives  very  accurate  results  (dependent  only  on  the  problem 
zoning)  when  checked  against  one -dimensional  problems  with  known  analyti¬ 
cal  solutions. 
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PARTICLE  VELOCITY  AT  SHOCK  FRONT  ON  CYLINDER  AXIS 

In  this  graph  (Figure  3),  the  shock  front  for  the  numerical  results  was  de¬ 
fined  +0  be  at  the  point  where  the  particle  velocity  is  a  maximum.  This 
graph  is  then  a  plot  of  maximum  particle  velocity  in  the  shock  versus  axial 
distance.  The  zero  of  axial  distance  has  been  adjusted  by  subtracting  5.  08 
cm  from  computational  values  in  order  to  compare  the  results  with  those  of 
NQLTR-bS-43.  ' 


Since  it  is  well  known  (see  Reference  8)  that  all  flow  quantities  overshoot 
their  true  values  when  the  Q-method  is  used,  these  results  should  be  some¬ 
what  high,  which  they  are.  A  more  accurate  method  of  the  determination  of 
the  particle  velocity  of  the  shock  front  would  involve  some  averaging.  Since 
such  techniques  are  somewhat  subjective  this  was  not  done,  Such  a  technique 
is,  however,  recommended  for  final  evaluation  of  this  data. 


Figure  3.  Particle  Velocity  at  Shock  Front  on  Cylinder  Axis 
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SHOCK  SHAPE 


The  shap°  of  the  shock  was  defined  here  by  the  position  of  maximum  Q.  These 

graphs  (Figures  4,  5  and  6)  are  thus  a  plot  of  position  of  ceils  of  maximum  Q 
in  the  shock.  No  zero  axial  distance  adjustment  was  made.  Figure  32  in 
Reference  4  shows  an  apparent  discontinuity  in  the  radius  of  curvature  of 
the  shock  wave  between  20  and  40  mm  in  the  Plexiglas.  Examination  of 
Figures  4,  5,  and  0  does  not  conclusively  demonstrate  such  behavior.  It  is 
the  opinion  of  these  authors  that  finer  zoning  could  add  much  to  the  resolution 
of  the  problem,  since  the  size  of  the  phenomena  in  question  is  comparable 
with  the  current  zone  size.  Again,  some  smoothing  of  the  data  would  be 
helpful.  This  was  not  done  because  of  the  subjective  nature  of  such  analysis. 


It  should  be  mentioned  that  an  abrupt  decrease  in  the  radius  of  curvature 
of  the  shock  wave  would  be  expected  when  the  outside  of  PMMA  material 
rarifies  through  20  kb.  .  At  this  point  the  sound  speed  would  decrease  dis- 
continuously,  and  sound  signals  would  thus  drive  the  shock  wave  with  less 
strength.  No  attempt  was  made  in  our  data  analysis  to  correlate  pressure 
levels  with  the  observed  discontinuity. 
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Figure  4.  Shock  Shape  at  8,  9,  10.  5  nsec 
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TIME  11.5h  SEC  TIME  12.5  mSEC  TIME  13.5  u  SEC 


Figure  5.  Shock  Shape  at  11.  5,  12.  5,  13.  5  psec 


TIME  14.5  MSEC  TIME  15.5  mSEC  TIME  16.5  „ SEC 


AXIAL  DJSTANrL  cm 

Figure  6.  Shock  Shape  at  14,  5,  15,  5,  16.  5  ^sec 


18 


NOLTR  69-219 


PARTICLE  VELOCITY 

Due  to  the  zoning  of  the  problem,  the  location  of  particles  of  interest  did 
not  fall  on  zoned  interfaces.  Thus,  particle  velocities  for  particles 
between  interfaces  had  to  be  determined  by  interpolation.  In  Figures  7 
through  21,  the  particle  velocities  at  the  interface  on  either  side  of  the 
location  of  interest  are  plotted.  A  linear  interpolation  is  made  between  these 
points  to  determine  the  particle  velocity  at  the  location  in  question.  No  zero 
time  adjustment  was  made.  No  attempt  was  made  to  correlate  particle  trajec¬ 
tories  with  elasto -plastic  wave  motion. 


Figure  7.  Particle  Velocity  versus  Time  (initial  axial 
location  -  0.  5  cm) 
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Figure  8.  Particle  Velocity  versus  Time  (initial  axial 
location  -  1.  0  cm) 


Figure  9.  Particle  Velocity  versus  Time  (initial  axial 
location  -  1.  5  cm) 
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Figure  10.  Particle  Velocity  versus  Time  (initial  axial 
location  --2.  0  cm) 


Figure  11.  Particle  Velocity  versus  Time  (initial  axial 
location  -  2.  5  cm) 
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Figure  16.  Particle  Velocity  versus  Time  (initial  axial 
location  -  5.  0  cm) 
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Figure  18.  Particle  Velocity  versus  Time  (initial  axial 
location  -  7.  0  cm) 
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AXIAL  PRESSURE  DISTRIBUTION 

This  graph  (Figure  22 1  is  typical  of  the  axial  pressure  distribution  across 
the  shock  front  in  the  Plexiglas.  It  shows  the  shockwave  and  following 
rarefaction  wave.  Similar  plots  of  the  detonation  wave  in  the  tetryl  indicated 
similar  behavior.  In  that  case,  peak  pressure  values  did  not  reach  the 
Chapman -Jouquet  value  due  to  the  extremely  sharp  (infinite  derivative)  de¬ 
crease  of  the  flow  quantities  from  the  Chapman-Jouquet  values  in  the  spherical 
Taylor  wave,  as  discussed  in  Reference  9.  Since  peak  shock  values  tend  to 
smooth  out  in  time,  it  is  believed  that  the  shock  values  in  the  Plexiglas  may  be 
quite  accurate.  Figure  2  3  tends  to  confirm  this  belief. 


Figure  22. 


Axial  Pressure  Distribution  at  the  Shock  I-  rent  m  the  PMMA 
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AXIAL  PRESSURE  AT  SHOCK  FRONT 

In  This  graph  (Figure  23)  the  shock  front  was  defined  to  be  at  the  location 
where  Q  is  a  maximum,.  A  zero  social  distance  adjustment  was  made  by 
subtracting  5,  08  cm  from  each  of  the  numerical  computed  axial  distance. 
Peak  values  of  pressure  which  follow  the  shock  front  as  defined  above  were 
plotted  as  the  true  pressure  values  at  the  shonk  front.  It  is  the  opinion  of 
these  authors  that  simulations  with  finer  zoning  would  be  required  if  a 
conclusive  comparison  of  shock  pressures  were  to  be  made.  This  is  espe¬ 
cially  true  in  the  region  from  0  to  4  cm. 


Figure  23.  Axial  Pressure  at  Shock  Front 
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SECTION  5 

DETAILED  DISPLAY  OF  THE  FLOW  DYNAMICS 

A  detailed  display  of  the  flow  dynamics  was  prepared  in  the  form  of  Calcomp 
pen  and  ink  plots  of  the  flow  field  quantities.  Because  of  their  bulk,  these 
are  gathered  in  Appendix  A  as  Figures  A1  through  A21  (representing  the  21 
simulation  time  intervals  selected  for  plotting  in  the  range  from  0.  01  psec 
to  32.  51  psec. 

Figures  identified  by  the  letter  (a),  e.  g. ;  Figure  Al(a),  are  plots  of  the 
Lagrangian  grid.  The  letter  (b)  following  the  figure  number  designates  plots 
of  particle  velocity,  and  (c)  indicates  plots  of  the  direction  of  maximum 
principal  stress,  with  stress  being  measured  positive  in  tension.  These 
figures  indicate  regions  of  particular  stress  patterns,  and  should  be  compared 
qualitatively  with  the  framing  camera  pictures  of  material  failure  in  Refer¬ 
ence  4.  As  discussed  pr°viously,  the  material  was  allowed  to  flow  plastically 
at  the  yield  strength  after  yield.  For  materials  such  as  Plexiglas  which  fail 
due  to  brittle  fracture  it  may  have  been  more  appropriate  to  set  all  the 
deviatoric  stresses  equal  to  zero  for  the  remainder  of  the  calculation,  and 
set  the  pressure  equal  to  zero  until  the  material  compresses  again,  if  it  does 
compress  again.  This  would  not  substantially  affect  the  particle  velocity. 

Figures  with  the  letter  (d)  are  isobar  plots.  Some  indication  of  shock  curva¬ 
ture  can  be  obtained  from  these.  Since  rarefaction  wave  fronts  are  character¬ 
istic  surfaces  of  the  partial  differential  equations  of  motion,  (i.  e.  ,  they  are 
surfaces  along  which  the  space  and  time  derivatives  of  the  flow  field  have 
discontinuities),  isobar  plots  should  have  discontinuous  slopes  at  a  wave 
front.  The  computation  performed  did  not  have  sufficient  resolution  to  indi¬ 
cate  such  wave  fronts.  In  examining  the  isobars  of  Reference  9  such  wave 
fronts  are  visible.  Those  calculations  were  of  course  done  with  much  more 
resolution.  This  again  is  another  indication  that  finer  zoned  simulations 
would  be  desirable. 
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Finally,  figures  with  a  letter  (e)  indicate  directions  of  maximum  principal 
strain  rate  weighted  with  the  magnitude  of  the  local  artificial  viscosity. 
Such  plots  indicate  positions  of  shock  and  detonation  waves. 

All  graphs  are  labeled  with  the  number  of  computer  integration  cycles 
(NCYCLE)  and  the  simulation  time  (TIME)  measured  from  the  time  the 
charge  of  1<etryl  was  initiated. 

It  should  be  remarked  at  this  point  that  all  of  the  Calcomp  plots  give  visual 
aid  in  understanding  the  flow,  but  do  not  replace  the  exact  quantitative 
information  available  from  the  computer  printouts.  Indeed,  the  Calcomps 
should  be  used  in  conjunction  with  the  printouts  for  future  quantitative 
evaluation  of  the  simulation. 
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SECTION  6 

CONCLUDING  REMARKS 

Results  of  the  simulation  of  the  shock  loading  of  PMMA  have  been  presented. 
It  is  felt  that  more  subjective  analysis  of  the  data,  which  would  include  such 
things  as  smoothing  and  interpolation,  would  be  helpful.  But,  since  these 
analyses  are  subjective,  they  were  not  done  here.  Beyond  this  it  would  be 
extremely  helpful  to  simulate  the  configuration  wit’  finer  zoning.  For  ex¬ 
ample,  the  shock  loading  at  the  PMMA  interface  could  be  simulated  with 
very  good  accuracy  by  using  the  analytical  Taylor  wave  for  that  portion  of  the 
"detonation  head"  which  is  bounded  by  the  leading  edge  of  the  detonation,  the 
lateral  rarefactions,  and  the  rarefactions  from  the  rear.  This  part  of  the 
detonation  wave  would  be  put  into  the  code  as  an  initial  condition.  It  would 
give  very  accurate  information  for  that  part  of  the  flow  which  is  determined 
by  the  appropriate  domain  of  dependance. 

A  separate  simulation  of  the  total  event,  done  with  resolution  at  least 
double  that  used  here  would  then  determine  events  farther  removed  from 
the  interface.  This  would  probably  determine  these  events  as  accurately  as 
much  higher  resolution  simulations,  due  to  the  lack  of  memory  of  hydro- 
dynamic  events.  Finally  it  is  recommended  that  a  serious  study  be  con¬ 
ducted  on  the  fracture  of  PMMA.  Such  a  study  would  involve  the  simulation 
of  several  constitutive  models  for  PMMA. 
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APPENDIX  A 

CALCOMP  PLOTS  OF  THE  FLOW  FIELD 
QUANTITIES  AT  VARIOUS  TIMES 


The  Calcomp  plots  shown  in  Figures  A1  through  A21  are  identified 
in  the  following  manner: 

Shown  on  the  face  of  each  plot  are  the  number  of  com¬ 
puter  integration  cycles  (NCYCLE)  and  the  simulation 
time  (TIME)  measured  from  initiation  of  detonation  in 
the  tetryl. 

Parenthetical  letters  following  the  figure  number,  e.g.  , 
Figure  Al(a),  indicate  plots  of: 

(a)  Lagrangian  grid 

(b)  Particle  velocity 

(c)  Direction  of  maximum  principal  stress 

(d)  Isobars 

(e)  Direction  of  maximum  principal  strain  rate  weighted 
with  the  magnitude  of  the  local  artificial  viscosity 

All  figures  do  not  necessarily  include  the  lull  range  of  plots,  (a) 
through  (e). 
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Figure  AG(b) 
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Figure  A 6(c) 
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Figure  A6(d) 
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